#####################
#####################

###Replication code for 
###Quotas and Party Priorities:
###Direct and Indirect Effects of Quota Laws 

###Main Analysis

#####################
#####################

rm(list = ls())

##set working directory
setwd("C:/Users/acw64/Dropbox/Working papers/Manifestos")

##load full data set

load("data_for_matching.RData")
head(prq_data) 

blah<-prq_data
nrow(blah) ##837

## Analysis with matched data 
############################
############################

load("matched_data.RData") ##data is dat2
head(dat2)
nrow(dat2) ##282

#################
####Table 1
#################

fm.1<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.1) 
fm.2<-lm(per504 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.2)  
fm.3<-lm(logrile ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.3) 

##install.packages("stargazer")
library(stargazer)

stargazer(fm.1, fm.2, fm.3, title="Regression Results",
align=TRUE, dep.var.labels=c("Social Justice","Welfare State Expansion", "Left-Right Position"),
omit.stat=c("LL","ser","f"), no.space=TRUE)

#################

##cluster SE function##

clx <- 
function(fm, dfcw, cluster){
         library(sandwich);library(lmtest)
         M <- length(unique(cluster))   
         N <- length(cluster)           
         K <- fm$rank                        
         dfc <- (M/(M-1))*((N-1)/(N-K))  
         uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
         vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
	return(list(coeftest= coeftest(fm, vcovCL), vcov=vcovCL))
 }

###############
### use to change SEs etc manually in latex

C1<- clx(fm.1, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.1$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.1$coefficients/std.err))),5), 
       lower=fm.1$coefficients-1.96*std.err, 
       upper=fm.1$coefficients+1.96*std.err)
       
r.est

C2<- clx(fm.2, 1, dat2$date)
std.err2<-sqrt(diag(C2$vcov))
r.est2<-cbind(estimate=fm.2$coefficients, std.err2, 
       pvalue=round(2*(1-pnorm(abs(fm.2$coefficients/std.err2))),5), 
       lower=fm.2$coefficients-1.96*std.err2, 
       upper=fm.2$coefficients+1.96*std.err2)
       
r.est2

C3<- clx(fm.3, 1, dat2$date)
std.err3<-sqrt(diag(C3$vcov))
r.est3<-cbind(estimate=fm.3$coefficients, std.err3, 
       pvalue=round(2*(1-pnorm(abs(fm.3$coefficients/std.err3))),5), 
       lower=fm.3$coefficients-1.96*std.err3, 
       upper=fm.3$coefficients+1.96*std.err3)
       
r.est3

#################
####Table 2
#################

fm.15<-lm(lag1.pfem_new ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.15) 

C15<- clx(fm.15, 1, dat2$date)
std.err15<-sqrt(diag(C15$vcov))
r.est15<-cbind(estimate=fm.15$coefficients, std.err15, 
       pvalue=round(2*(1-pnorm(abs(fm.15$coefficients/std.err15))),5), 
       lower=fm.15$coefficients-1.96*std.err15, 
       upper=fm.15$coefficients+1.96*std.err15)
       
r.est15

fm.16<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.pfem_new + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.16) 

C16<- clx(fm.16, 1, dat2$date)
std.err16<-sqrt(diag(C16$vcov))
r.est16<-cbind(estimate=fm.16$coefficients, std.err16, 
       pvalue=round(2*(1-pnorm(abs(fm.16$coefficients/std.err16))),5), 
       lower=fm.16$coefficients-1.96*std.err16, 
       upper=fm.16$coefficients+1.96*std.err16)
       
r.est16

stargazer(fm.15, fm.16, title="Regression Results",
align=TRUE, dep.var.labels=c("Women in Party (t-1)","Social Justice"),
omit.stat=c("LL","ser","f"), no.space=TRUE)

###again change SE with clustered from clx function above

#################
####Figure 2
#################

library(sandwich)
library(lmtest)

fm.17<-lm(per503 ~ factor(party) + factor(year) + quota_ey3_v2 + quota_ey2_v2 + quota_ey1_v2 + quota_lead1_v2 + quota_lead2_v2 + quota_lead3_v2  
+ lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.17) ##
C1<- clx(fm.17, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.17$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.17$coefficients/std.err))),5), 
       lower=fm.17$coefficients-1.96*std.err, 
       upper=fm.17$coefficients+1.96*std.err)
       
r.est

###model coefs, lb and ub saved in csv file, "data_leads_lags_per503_v3.csv"

##########Make the figure 

mydata4 = read.csv("data_leads_lags_per503_v3.csv")
head(mydata4)


library(ggplot2)
library(gridExtra)

bp<-ggplot(mydata4, aes( x = factor(model2), y = coef)) +
  geom_point(size = 3)+ geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymax = ub, ymin = lb), size=1, width=.5) + 
  ylab("Coefficient (95% CIs)") + xlab("") +  theme_bw()  + theme(legend.position="none") + theme(text = element_text(size=20))



bp2<- bp + scale_x_discrete(labels=c("1" = "3 elections prior",
                              "2" = "2 elections prior", "3" = "1 election prior", "4" = "Election quota first implemented",
"5" = "1 election after", "6" = "2 elections after"))  + theme(
axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'))
bp2

#####################
#####################

## Appendix 
############################
############################

##quota shock variable is in full data set; merge to incude in matched data 
#covariates<- c("party", "date", "quota_shock", "lag1.quota_shock")
#X3<-dat[,covariates]
#tail(X3)

#test5<- merge(dat2, X3, by=c("party", "date"), all.x=TRUE) #no NAs
#head(test5)
#nrow(test5)
#nrow(dat2) 
#dat2<-test5

##and same for full 

#test6<- merge(blah, X3, by=c("party", "date"), all.x=TRUE) #no NAs
#head(test6)
#nrow(test6)
#nrow(blah) 
#blah<-test6

############################
###Table B1: Summary statistics, matched data 
library(pastecs)
stat.desc(dat2)

############################
###Table B2: Summary statistics, full data 
stat.desc(blah)

############################
###Table C1: Moderated Effects of Quota Laws on Policy Positions


fm.17<-lm(per503 ~ factor(party) + factor(year) +  lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap) + lag1.quota + lag1.quota:right, data=dat2)
summary(fm.17) 


###Wald Test - quota and interaction term 
library(aod)

l <- cbind(0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,-1)

wald.test(b=coef(fm.17), Sigma=vcov(fm.17), L=l) ##p = 0.000

C17<- clx(fm.17, 1, dat2$date)
std.err17<-sqrt(diag(C17$vcov))
r.est17<-cbind(estimate=fm.17$coefficients, std.err17, 
       pvalue=round(2*(1-pnorm(abs(fm.17$coefficients/std.err17))),5), 
       lower=fm.17$coefficients-1.96*std.err17, 
       upper=fm.17$coefficients+1.96*std.err17)
       
r.est17


###NOTE: interaction also not sig for other DVs: 
fm.17b<-lm(per504 ~ factor(party) + factor(year) + lag1.quota + lag1.quota:right + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.17b)
fm.17c<-lm(rile ~ factor(party) + factor(year) + lag1.quota + lag1.quota:right + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.17c) ##also tried logrile


############################
##Figure C1: Predicted Change in Social Justice as a Function of Quota Law and Party Type


###plotting the interaction 
library(interplot)
library(ggplot2)
library("ggthemes")
library("scales")

p1<-interplot(m = fm.17, var1 = "lag1.quota", var2 = "right")+ 
 ylab("Estimated Coefficient for Social Justice") + geom_point(size = 3)  +
  geom_errorbar(aes(ymax = ub, ymin = lb), size=1, width=0.3) + 
 theme_classic() + theme(legend.position="none") + theme(text = element_text(size=18)) + 
  scale_fill_brewer(palette = "Set1")

p2<-p1 + scale_x_discrete(limit = c(0,1),
                     labels = c("Left Party","Right Party"))


p2

##############################
### Check for interaction between party quota and quota law


fm.313<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.quota:lag1.partyquota  + lag1.fplabfo  + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.313) 


#################

C1<- clx(fm.313, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.313$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.313$coefficients/std.err))),5), 
       lower=fm.313$coefficients-1.96*std.err, 
       upper=fm.313$coefficients+1.96*std.err)
       
r.est


###plotting the interaction (Figure C2)
library(interplot)
library(ggplot2)
library("ggthemes")
library("scales")

p1<-interplot(m = fm.313, var1 = "lag1.quota", var2 = "lag1.partyquota")+ 
 ylab("Estimated Coefficient for Social Justice") + geom_point(size = 3)  +
  geom_errorbar(aes(ymax = ub, ymin = lb), size=1, width=0.3) + 
 theme_classic() + theme(legend.position="none") + theme(text = element_text(size=18)) + 
  scale_fill_brewer(palette = "Set1")

p2<-p1 + scale_x_discrete(limit = c(0,1),
                     labels = c("No Party Quota","Party Quota"))
p2

############################
##Table D1: Regression Results, Full Data Set


fm.8<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=blah)
summary(fm.8) ##results the same with no controls 
fm.9<-lm(per504 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=blah)
summary(fm.9) ##
fm.10<-lm(logrile ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=blah)
summary(fm.10)

stargazer(fm.8, fm.9, fm.10, title="Regression Results",
align=TRUE, dep.var.labels=c("Social Justice","Welfare State Expansion", "Left-Right Position"),
omit.stat=c("LL","ser","f"), no.space=TRUE)


C8<- clx(fm.8, 1, blah$date)
std.err8<-sqrt(diag(C8$vcov))
r.est8<-cbind(estimate=fm.8$coefficients, std.err8, 
       pvalue=round(2*(1-pnorm(abs(fm.8$coefficients/std.err8))),5), 
       lower=fm.8$coefficients-1.96*std.err8, 
       upper=fm.8$coefficients+1.96*std.err8)
       
r.est8

C9<- clx(fm.9, 1, blah$date)
std.err9<-sqrt(diag(C9$vcov))
r.est9<-cbind(estimate=fm.9$coefficients, std.err9, 
       pvalue=round(2*(1-pnorm(abs(fm.9$coefficients/std.err9))),5), 
       lower=fm.9$coefficients-1.96*std.err9, 
       upper=fm.9$coefficients+1.96*std.err9)
       
r.est9

C10<- clx(fm.10, 1, blah$date)
std.err10<-sqrt(diag(C10$vcov))
r.est10<-cbind(estimate=fm.10$coefficients, std.err9, 
       pvalue=round(2*(1-pnorm(abs(fm.10$coefficients/std.err10))),5), 
       lower=fm.10$coefficients-1.96*std.err10, 
       upper=fm.10$coefficients+1.96*std.err10)
       
r.est10

############################
##Table D2: Regression Results, Using �Quota Impact� Variable

fm.18b<-lm(per503 ~ factor(party) + factor(year) + lag1.quota_shock + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.18b) ##

C16<- clx(fm.18b, 1, dat2$date)
std.err16<-sqrt(diag(C16$vcov))
r.est16<-cbind(estimate=fm.18b$coefficients, std.err16, 
       pvalue=round(2*(1-pnorm(abs(fm.18b$coefficients/std.err16))),5), 
       lower=fm.18b$coefficients-1.96*std.err16, 
       upper=fm.18b$coefficients+1.96*std.err16)
       
r.est16

fm.2b<-lm(per504 ~ factor(party) + factor(year) + lag1.quota_shock + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.2b) ##results are the same when no controls are included 

C16<- clx(fm.2b, 1, dat2$date)
std.err16<-sqrt(diag(C16$vcov))
r.est16<-cbind(estimate=fm.2b$coefficients, std.err16, 
       pvalue=round(2*(1-pnorm(abs(fm.2b$coefficients/std.err16))),5), 
       lower=fm.2b$coefficients-1.96*std.err16, 
       upper=fm.2b$coefficients+1.96*std.err16)
       
r.est16

fm.3b<-lm(logrile ~ factor(party) + factor(year) + lag1.quota_shock + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.3b) ##p = 0.09 with no controls

C16<- clx(fm.3b, 1, dat2$date)
std.err16<-sqrt(diag(C16$vcov))
r.est16<-cbind(estimate=fm.3b$coefficients, std.err16, 
       pvalue=round(2*(1-pnorm(abs(fm.3b$coefficients/std.err16))),5), 
       lower=fm.3b$coefficients-1.96*std.err16, 
       upper=fm.3b$coefficients+1.96*std.err16)
       
r.est16


############################
##Table D3: Regression Results, with LDV (no FE)

fm.5<-lm(per503 ~ lag1.per503 + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.5) ##with LDV alone (and trend), sig, also when including party FE (with or without LDV), and party and year FE
##  1.91363 , p= 0.0226 *

C1<- clx(fm.5, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.5$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.5$coefficients/std.err))),5), 
       lower=fm.5$coefficients-1.96*std.err, 
       upper=fm.5$coefficients+1.96*std.err)
       
r.est

fm.6<-lm(per504 ~ lag1.per504 + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.6) 

C1<- clx(fm.6, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.6$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.6$coefficients/std.err))),5), 
       lower=fm.6$coefficients-1.96*std.err, 
       upper=fm.6$coefficients+1.96*std.err)
       
r.est

fm.7<-lm(logrile ~ lag1.logrile + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.7)

C1<- clx(fm.7, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.7$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.7$coefficients/std.err))),5), 
       lower=fm.7$coefficients-1.96*std.err, 
       upper=fm.7$coefficients+1.96*std.err)
       
r.est


############################
##Table D4: Regression Results, No Control Variables (FE Retained)

fm.1a<-lm(per503 ~ factor(party) + factor(year) + lag1.quota, data=dat2)
summary(fm.1a) ##results are the same when no controls are included 
fm.2a<-lm(per504 ~ factor(party) + factor(year) + lag1.quota, data=dat2)
summary(fm.2a) ##results are the same when no controls are included 
fm.3a<-lm(logrile ~ factor(party) + factor(year) + lag1.quota, data=dat2)
summary(fm.3a) ##p = 0.09 with no controls


C1<- clx(fm.1a, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.1a$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.1a$coefficients/std.err))),5), 
       lower=fm.1a$coefficients-1.96*std.err, 
       upper=fm.1a$coefficients+1.96*std.err)
       
r.est

C2<- clx(fm.2a, 1, dat2$date)
std.err2<-sqrt(diag(C2$vcov))
r.est2<-cbind(estimate=fm.2a$coefficients, std.err2, 
       pvalue=round(2*(1-pnorm(abs(fm.2a$coefficients/std.err2))),5), 
       lower=fm.2a$coefficients-1.96*std.err2, 
       upper=fm.2a$coefficients+1.96*std.err2)
       
r.est2

C3<- clx(fm.3a, 1, dat2$date)
std.err3<-sqrt(diag(C3$vcov))
r.est3<-cbind(estimate=fm.3a$coefficients, std.err3, 
       pvalue=round(2*(1-pnorm(abs(fm.3a$coefficients/std.err3))),5), 
       lower=fm.3a$coefficients-1.96*std.err3, 
       upper=fm.3a$coefficients+1.96*std.err3)
       
r.est3



############################
##Table D5: Regression Results, Lagging Control Variables by 2 Election-Years


covariates<- c("party", "date", "lag2.partyquota", "lag2.fplabfo", "lag2.pervote", "lag2.rgdpecap", "lag2.effpar_leg")
X4<-dat[,covariates] 
head(X4)

test6<- merge(dat2, X4, by=c("party", "date"), all.x=TRUE) #no NAs
head(test6)
nrow(test6)
nrow(dat2) 

dat2<-test6
testing<-na.omit(dat2)

fm.32<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag2.partyquota + lag2.fplabfo + lag2.pervote + lag2.effpar_leg + log(lag2.rgdpecap), data=testing)
summary(fm.32) ##results similar, coef = 5.65 **

C1<- clx(fm.32, 1, testing$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.32$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.32$coefficients/std.err))),5), 
       lower=fm.32$coefficients-1.96*std.err, 
       upper=fm.32$coefficients+1.96*std.err)
       
r.est

fm.33<-lm(per504 ~ factor(party) + factor(year) + lag1.quota + lag2.partyquota  + lag2.fplabfo + lag2.pervote + lag2.effpar_leg + log(lag2.rgdpecap), data=testing)
summary(fm.33) ##still no effect

C1<- clx(fm.33, 1, testing$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.33$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.33$coefficients/std.err))),5), 
       lower=fm.33$coefficients-1.96*std.err, 
       upper=fm.33$coefficients+1.96*std.err)
       
r.est

fm.34<-lm(logrile ~ factor(party) + factor(year) + lag1.quota + lag2.partyquota + lag2.fplabfo + lag2.pervote + lag2.effpar_leg + log(lag2.rgdpecap), data=testing)
summary(fm.34) ##still no effect 

C1<- clx(fm.34, 1, testing$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.34$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.34$coefficients/std.err))),5), 
       lower=fm.34$coefficients-1.96*std.err, 
       upper=fm.34$coefficients+1.96*std.err)
       
r.est


############################
##Table D6: Regression Results, Excluding Parties that Proposed a Quota

##now try taking out parties that supported a quota law. 
tb <- subset(dat2, party!=21321 & party!=33320 & party!=31320 & party!=35311 & party!=21521 )
nrow(tb) ##242
table(tb$quota, tb$left)

#dat2$partyname[dat2$party==21221]
#table(dat2$partyname[dat2$countryname=="France"], dat2$party[dat2$countryname=="France"])

##33320 Spanish Socialist Workers’ Party
##31320 French Socialist Party ##NOTE that the RPR which also supported quota (31625) is not in the data
# because it only existed in the 1997 -- 2002 period. 
##35311 Portugal Socialist Party
##21521 Belgium Christian Democratic and Flemish
##21321 Flemish Socialist Party 

fm.21<-lm(per503 ~  factor(party) + factor(year) + lag1.quota + 
 lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=tb)
summary(fm.21) ##results slightly increase in size 

C21<- clx(fm.21, 1, tb$date)
std.err21<-sqrt(diag(C21$vcov))
r.est21<-cbind(estimate=fm.21$coefficients, std.err21, 
       pvalue=round(2*(1-pnorm(abs(fm.21$coefficients/std.err21))),5), 
       lower=fm.21$coefficients-1.96*std.err21, 
       upper=fm.21$coefficients+1.96*std.err21)
       
r.est21 ##5.916,p = 0.000
## NOTE that there is also no change to other models

fm.22<-lm(per504 ~  factor(party) + factor(year) + lag1.quota + 
 lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=tb)
summary(fm.22)

C22<- clx(fm.22, 1, tb$date)
std.err22<-sqrt(diag(C22$vcov))
r.est22<-cbind(estimate=fm.22$coefficients, std.err22, 
       pvalue=round(2*(1-pnorm(abs(fm.22$coefficients/std.err22))),5), 
       lower=fm.22$coefficients-1.96*std.err22, 
       upper=fm.22$coefficients+1.96*std.err22)
       
r.est22 ##

fm.23<-lm(logrile ~  factor(party) + factor(year) + lag1.quota + 
 lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=tb)
summary(fm.23)

C23<- clx(fm.23, 1, tb$date)
std.err23<-sqrt(diag(C23$vcov))
r.est23<-cbind(estimate=fm.23$coefficients, std.err23, 
       pvalue=round(2*(1-pnorm(abs(fm.23$coefficients/std.err23))),5), 
       lower=fm.23$coefficients-1.96*std.err23, 
       upper=fm.23$coefficients+1.96*std.err23)
       
r.est23 ##

############################
##Table D7: Regression Results, Unlagged variables

##merge some of these from original data set 

covariates<- c("party", "date", "fplabfo", "effpar_leg", "rgdpecap")
X3<-dat[,covariates]
tail(X3)

test5<- merge(dat2, X3, by=c("party", "date"), all.x=TRUE) #no NAs
head(test5)
nrow(test5)
nrow(dat2) 
dat2<-test5


fm.15<-lm(pfem_new ~ factor(party) + factor(year) + quota + partyquota + + fplabfo + pervote + effpar_leg + log(rgdpecap)  , data=dat2)
summary(fm.15) ##it doesn't matter (interesting to note that party quota is positive)

C23<- clx(fm.15, 1, dat2$date)
std.err23<-sqrt(diag(C23$vcov))
r.est23<-cbind(estimate=fm.15$coefficients, std.err23, 
       pvalue=round(2*(1-pnorm(abs(fm.15$coefficients/std.err23))),5), 
       lower=fm.15$coefficients-1.96*std.err23, 
       upper=fm.15$coefficients+1.96*std.err23)
       
r.est23 ##quotas are not linked to increase in women's rep at party level

fm.165<-lm(per503 ~ factor(party) + factor(year) + quota + partyquota + + fplabfo + pervote + effpar_leg + log(rgdpecap)  , data=dat2)
summary(fm.165) 


C24<- clx(fm.165, 1, dat2$date)
std.err23<-sqrt(diag(C24$vcov))
r.est24<-cbind(estimate=fm.165$coefficients, std.err23, 
       pvalue=round(2*(1-pnorm(abs(fm.165$coefficients/std.err23))),5), 
       lower=fm.165$coefficients-1.96*std.err23, 
       upper=fm.165$coefficients+1.96*std.err23)
       
r.est24

stargazer(fm.15, fm.165, title="Regression Results",
align=TRUE, dep.var.labels=c("Women in Party","Social Justice"),
omit.stat=c("LL","ser","f"), no.space=TRUE)

############################
##Table D8: Tests of common trends assumption

fm.17<-lm(per503 ~ factor(party) + factor(year) + quota_ey3_v2 + quota_ey2_v2 + quota_ey1_v2 + quota_lead1_v2 + quota_lead2_v2 + quota_lead3_v2  
+ lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.17) ##
C1<- clx(fm.17, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.17$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.17$coefficients/std.err))),5), 
       lower=fm.17$coefficients-1.96*std.err, 
       upper=fm.17$coefficients+1.96*std.err)
       
r.est

stargazer(fm.17, title="Regression Results",
align=TRUE, dep.var.labels=c("Social Justice"),
omit.stat=c("LL","ser","f"), no.space=TRUE)

#################
##Party specific time trends
fm.111<-lm(per503 ~ factor(party)*year + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.111) 

C3<- clx(fm.111, 1, dat2$date)
std.err3<-sqrt(diag(C3$vcov))
r.est3<-cbind(estimate=fm.111$coefficients, std.err3, 
       pvalue=round(2*(1-pnorm(abs(fm.111$coefficients/std.err3))),5), 
       lower=fm.111$coefficients-1.96*std.err3, 
       upper=fm.111$coefficients+1.96*std.err3)
       
r.est3 ##results now not sig

############################
####Other tests cited in text 

############################
##robustness check: the other left-right scale 
############################

fm.11<-lm(rile ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.11) 


############################
##Robustness check: Full data set moderation (women in party) analysis
############################

fm.16<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.pfem_new + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=blah)
summary(fm.16) ##no change in results 


##########################
## Robustness checks: taking out quota countries one by one 
##########################

coun <- subset(dat2, countryname!="Spain")
nrow(coun) ##
table(coun$quota)
##repeat for France, Belgium, Portugal 

## Run baseline model again
fm.28<-lm(per503 ~ factor(party) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=coun)
summary(fm.28) 
## when remove belgium, coef drops to 3.5 and only sig at .1 level (p=0.08)

C28<- clx(fm.28, 1, coun$date)
std.err28<-sqrt(diag(C28$vcov))
r.est28<-cbind(estimate=fm.28$coefficients, std.err28, 
       pvalue=round(2*(1-pnorm(abs(fm.28$coefficients/std.err28))),5), 
       lower=fm.28$coefficients-1.96*std.err28, 
       upper=fm.28$coefficients+1.96*std.err28)
       
r.est28

############################
##Robustness check: Heterogeneous effects by country 


fm.15<-lm(pfem_new ~ factor(countryname) + factor(year) + lag1.quota + lag1.quota:factor(countryname) +
 lag1.partyquota + lag1.fplabfo + lag1.pervote + lag1.effpar_leg + log(lag1.rgdpecap), data=dat2)
summary(fm.15) ##it doesn't matter (interesting to note that party quota is positive)

C23<- clx(fm.15, 1, dat2$date)
std.err23<-sqrt(diag(C23$vcov))
r.est23<-cbind(estimate=fm.15$coefficients, std.err23, 
       pvalue=round(2*(1-pnorm(abs(fm.15$coefficients/std.err23))),5), 
       lower=fm.15$coefficients-1.96*std.err23, 
       upper=fm.15$coefficients+1.96*std.err23)
       
r.est23 ##nothing significant


##############################
### Country FE instead of Party FE 

fm.31<-lm(per503 ~ factor(countryname) + factor(party) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + log(lag1.rgdpecap), data=dat2)
summary(fm.31)  ##many parties drop out of the analysis when both party and country FE are included 


fm.31<-lm(per503 ~ factor(countryname) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote + log(lag1.rgdpecap), data=dat2)
summary(fm.31) 
fm.32<-lm(per504 ~ factor(countryname) + factor(year) + lag1.quota + lag1.partyquota  + lag1.fplabfo + lag1.pervote  + log(lag1.rgdpecap), data=dat2)
summary(fm.32)  
fm.33<-lm(logrile ~ factor(countryname) + factor(year) + lag1.quota + lag1.partyquota + lag1.fplabfo + lag1.pervote + log(lag1.rgdpecap), data=dat2)
summary(fm.33) 

library(stargazer)

stargazer(fm.31, fm.32, fm.33, title="Regression Results",
align=TRUE, dep.var.labels=c("Social Justice","Welfare State Expansion", "Left-Right Position"),
omit.stat=c("LL","ser","f"), no.space=TRUE)

#################

C1<- clx(fm.31, 1, dat2$date)
std.err<-sqrt(diag(C1$vcov))
r.est<-cbind(estimate=fm.31$coefficients, std.err, 
       pvalue=round(2*(1-pnorm(abs(fm.31$coefficients/std.err))),5), 
       lower=fm.31$coefficients-1.96*std.err, 
       upper=fm.31$coefficients+1.96*std.err)
       
r.est

C2<- clx(fm.32, 1, dat2$date)
std.err2<-sqrt(diag(C2$vcov))
r.est2<-cbind(estimate=fm.32$coefficients, std.err2, 
       pvalue=round(2*(1-pnorm(abs(fm.32$coefficients/std.err2))),5), 
       lower=fm.32$coefficients-1.96*std.err2, 
       upper=fm.32$coefficients+1.96*std.err2)
       
r.est2

C3<- clx(fm.33, 1, dat2$date)
std.err3<-sqrt(diag(C3$vcov))
r.est3<-cbind(estimate=fm.33$coefficients, std.err3, 
       pvalue=round(2*(1-pnorm(abs(fm.33$coefficients/std.err3))),5), 
       lower=fm.33$coefficients-1.96*std.err3, 
       upper=fm.33$coefficients+1.96*std.err3)
       
r.est3




